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The overall structure of the transition state and intermediate ensembles experimentally observed 
for Dihydrofolate Reductase and Interleukin-1(3 can be obtained utilizing simplified models which 
have almost no energetic frustration. The predictive power of these models suggest that, even for 
these very large proteins with completely different folding mechanisms and functions, real protein 
sequences are sufficiently well designed and much of the structural heterogeneity observed in the 
intermediates and the transition state ensembles is determined by topological effects. 
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I. INTRODUCTION 

Explaining how proteins self-assemble into well defined 
structures is a longstanding challenge. Energy land- 
scape theory and the funnel concept (|l|-0) have pro- 
vided the theoretical framework necessary for improving 
our understanding of this problem — efficient folding se- 
quences minimize frustration. Frustration may arise from 
the inability to satisfy all native interactions and from 
strong non-native contacts which can create conforma- 
tional traps. The difficulty of minimizing energetic frus- 
tration by sequence design, however, is also dependent 
on the choice of folding motif. Some folding motifs are 
easier to design than others (||,^) , suggesting the possibil- 
ity that evolution not only selected sequences with suffi- 
ciently small energetic frustration but also selected more 
easily designable native structures. To address this dif- 
ference in foldability, we have introduced the concept of 
"topological frustration" (10-13) — even when sequences 
have been designed with minimal energetic frustration, 
variations in the degree of nativeness of contacts in the 
transition state ensemble (TSE) are observed because of 
asymmetries imposed by the chosen final structure. 

Recent theoretical and experimental evidences sug- 
gest that proteins, especially small fast folding (sub- 
millisecond) proteins, have sequences with a sufRciently 
reduced level of energetic frustration that the global char- 
acteristics of the observed heterogeneity observed in the 
TSE are strongly influenced by the native state topol- 
ogy. We have shown (^3|) that the overall structure of 
the TSE for Chymotrypsin Inhibitor 2 (CI2) and for the 
SH3 domain of the src tyro sine-protein kinease can be 
obtained by using simplified models constructed by us- 
ing sequences that have almost no energetic frustration 
(Go-like potentials). These models drastically reduce the 
energetic frustration and energetic heterogeneity for na- 
tive contacts, leaving the topology as the primary source 
of the residual frustration. Topological effects, however, 
go beyond affecting the structure of the TSE. The over- 
all structure of the populated intermediate state ensem- 
bles during the folding of proteins such as Barnase, Ri- 
bonuclease H and Che Y have also been successfully de- 
termined using a similar model ([T3|). It is interesting 
to notice that although these model, since they consider 
totally unfrustrated sequences, may not reproduce the 
precise energetics of the real proteins, such as the value 
of the barrier heights and the stability of the intermedi- 
ates, they are able to determine the general structure of 
these ensembles. Therefore, the fact that these almost 
energetically unfrustrated models reproduce most of the 
major features of the TSE of these proteins indicate that 
real protein sequences are sufficiently well designed (i.e. 
with reduced energetic frustration) that much of the het- 
erogeneity observed in the TSE's and intermediates have 
a strong topological dependence. 

Do these conclusions hold to larger and slower fold- 
ing proteins with a more complex folding kinetics than 



two-state folders as CI2 and SH3? The success obtained 
with Barnase, Ribonuclease H and C/ie F intermediates 
already provides some encouragement — topology ap- 
pears to be important in determining on-pathway fold- 
ing intermediates. In this paper this approach is ex- 
tended to a pair of larger proteins: Dihydrofolate Re- 
ductase (DHFR) and Interleukin-ip (IL-1/3). The syn- 
optic analysis of these two proteins is particularly inter- 
esting because they have a comparable size (slightly over 
150 amino-acids), but different native structures, folding 
mechanisms and functions: DHFR is a two-domain a/P 
enzyme that maintains pools of tetrahydrofolate used in 
nucleotide metabolism while IL-1/3 is a single domain all 
/3 cytokine with no catalytic activity on its own but elicits 
a biological response by binding to its receptor. 



II. NUMERICAL PROCEDURES 

The energetically unfrustrated model of DHFR and 
IL-1/3 are constructed by using a Go-like Hamiltonian 
(|l^,^5|). A Go-like potential takes into account only na- 
tive interactions, and each of these interactions enters in 
the energy balance with the same weight. Residues in 
the proteins are represented as single beads centered in 
their C-a positions. Adjacent beads are strung together 
into a polymer chain by means of bond and angle inter- 
actions, while the geometry of the native state is encoded 
in the dihedral angle potential and a non-local bead-bead 
potential. 

A detailed description of this energy function can be 
found elsewhere (|l^). The local (torsion) and non-local 
terms have been adjust so that the stabilization energy 
residing in the tertiary contacts is approximately twice as 
large as the torsional contribution. This balance among 
the energy terms is optimal for the folding of our Go-like 
protein models (^. Solvent mediation and side chain ef- 
fects are already included in these effective energy func- 
tions. Therefore, entropy changes are associated to the 
configurational entropy of the chain. The native con- 
tact map of a protein is derived with the CSU software. 



based upon the approach developed in ref. (16|). Native 
contacts between pairs of residues {i,j) with j < i -\- A 
are discarded from the native map as any three and four 
subsequent residues are already interacting in the angle 
and dihedral terms. A contact between two residues (z, j) 
is considered formed if the distance between the C^s is 
shorter than 7 times their native distance CTy . It has been 
shown (|ri|) that the results are not strongly dependent 
on the choice made for the cut-off distance 7. In this 
work we used 7 — 1.2. 

For both (DHFR and IL-1/3) protein models, folding 
and unfolding simulations have been performed at sev- 
eral temperatures around the folding temperature. The 
results from the different simulations have been com- 
bined using the WHAM algorithm ([l7|). Several very 
different initial unfolded structures for the folding simu- 
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lations have been selected and they have been obtained 
from high temperature unfolding simulations. In order to 
have appropriate statistics, we made sure that for every 
transition state ensemble or intermediate, we have sam- 
pled about 500 uncorrelated conformations (thermally 
weighted). For smaller proteins such as SH3 and CI2 
(that have about 1/3 of the tertiary contacts of DHFR 
and IL-1/3) we have determined that about 200 uncorre- 
lated conformations in the transition state ensemble are 
necessary to have an error on the estimates of contact 
probabilities (or <I> values) of ± 0.05 (^. 

III. COMPARING SIMULATIONS AND 
EXPERIMENTS FOR DIHYDROFOLATE 
REDUCTASE AND INTERLEUKIN-1/? 

Dihydrofolate Reductase and Interleukin-1(3 not only 
have dissimilar native fold^ but also the nature of the 
intermediate states populated during the folding event is 
remarkably different. To explore the connection between 
the protein topology and the nature of the intermedi- 
ates, we used an energetically minimally frustrated 
model for these two proteins, with a potential energy 
function defined by considering only the native local and 
non-local interactions as being attractive (see Numerical 
Procedures, for details). This is a very simplified poten- 
tial that retains only information about the native fold 
— energetic frustration is almost fully removed. Notice 
that although the real amino-acid sequence is not in- 
cluded in this model, the chosen potential is like a "per- 
fect" sequence for the target structure, without the en- 
ergetic frustration of real sequences (since this potential 
includes attractive native tertiary contacts, it implicitly 
incorporates hydrophobic interactions). Therefore, this 
model provide us with the perfect computational tool to 
investigate how much of the structural heterogeneity ob- 
served during folding mechanism could be inferred from 
the knowledge of the native structure alone without con- 
tributions from energetic frustration. 

Since early work suggests that proteins (at least small 
fast folding proteins) have sufficiently reduced energetic 
frustration, they have a funnel-like energy landscape with 
a solvent-averaged potential strongly correlated with the 
degree of nativeness (but with some roughness due to the 
residual frustration). In this situation, the folding dy- 
namics can be described as the diffusion of an ensemble 
of protein configurations over a low dimensional free en- 
ergy surface — defined in terms of the reaction coordinate 



'^The 162 residues of DHFR arrange themselves in 8 (3- 
strands and 4 Q-helices, grouped together in the folded state 
in as detailed in Fig.|| (d), while IL-1/3 is a 153 residues, all-/? 
protein, composed by 12 /3 strands packed together as shown 
in Figg (c)-(d). 



Q, where Q represents the fraction of the native contact 
formed in a conformation (Q = at the fully unfolded 
state and Q = 1 at the folded state) (|l§^,^. The 
ensemble of intermediates observed in this free energy 
profile are expected to mimic the real kinetic intermedi- 
ates. 

Fig. |l| shows a comparison between the folding mech- 
anism obtained from our simulations for the minimally- 
frustrated analogue of DHFR (panels (a) and (c)) and 
IL-1/? (panels (b) and (d)). The different nature of the 
folding intermediates of the two proteins and their na- 
tive ensembles emerging from these data is in substantial 
agreement with the experimental observations, with the 
adenine binding domain of DHFR being folded in the 
main intermediate in the simulation and the central (3 
strands of IL-1/3 being formed early in this single domain 
protein. The absolute values of the free energy barri- 
ers resulting from simulations may not necessarily agree 
with the experimental ones because we are dealing with 
unfrustrated designed sequences. Thus, quantitative pre- 
dictions that depend on barrier heights and stability of 
the intermediate ensembles (e.g. folding time, rate de- 
termining barriers and lifetime of intermediates) are not 
possible for this kind of models. However we show that 
topology is sufficient to correctly detect the positions of 
the transition state and intermediate states. A more de- 
tailed description follows. 

A. Dihydrofolate Reductase 

The folding process emerging from the dynamics of 
the Go-like analogue of DHFR (as summarized in Fig. 
|l| (a) and (c)) is interestingly peculiar and consistent with 
the experimentally proposed folding mechanism ( p^ (see 
Fig. |^(d)). Refolding initiates by a barrier less collapse to 
a quasi-stable species (Q=0.2) which corresponds to the 
formation of a burst-phase intermediate, I bp, with little 
stability but some protection from iJ-exchange across the 
central (3 sheet (pO|). This initial collapse is followed by 
production of the main intermediate Ihf (Highly Flu- 
orescent), which is described in the mechanism of Fig. 
^ (d) as the collection of intermediates Ii — I4. Ii — I4 
are structurally similar to each other, but differentiated 
experimentally by the rate at which they proceed to- 
wards the native protein. Finally, after the overcom- 
ing of a second barrier, the protein visits an ensemble 
of native structures with different energies. The experi- 
mentally determined folding mechanism of DHFR shows 
transient kinetic control in the formation of native con- 
formers (iV4 dominant). This is later overridden by ther- 
modynamic considerations {N2 dominant) at final equi- 
librium (p^. This latter finding is consistent with the 
nature of the folding ensemble determined by the simu- 
lations. As shown in Fig. |^ (b) a set of structures close 
to to the native state (Q around 0.7-0.8) is transiently 
populated beside the fully folded state (Q = 1). Since 
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the main intermediate Ihf has been recently character- 
ized by experimental studies (|l9|-^^, we take our anal- 
ysis a step farther by comparing the average structure 
of the Ih f ensemble from our simulations to the one ex- 
perimentally determined. For this purpose we compute 
the formation probability Qij{Q) for each native DHFR 
contact -involving residues at different stages of 

the folding process by averaging the number of times the 
contact occurs over the set of structures existent in a se- 
lected range of Q- As detailed in Fig. ||, the central result 
from this analysis is that the main intermediate Ihf is 
characterized by a largely different degree of formation 
in different parts of the protein: domain 1 (i.e. interac- 
tions among strands 2-5 and helices 2-3) appears to be 
formed with probability greater than 0.7 while domain 2 
(i.e. interactions among strands 6-8, helix 1 and helix 4) 
is almost non existent. 

The formation of domain 1 and domain 2 during the 
folding event is more closely understood from Fig. |^ (pan- 
els (a) and (c)), where the RMS distance of the parts 
of the protein constituting each domain from the corre- 
sponding native structures is shown for a typical folding 
simulation. Indeed the two domains fold in a noticeably 
different way: in the stable intermediate Ihf, domain 
1 is closer than 5 A (RMS) to that found in the native 
structure while domain 2 is highly variable (RMS dis- 
tance greater than 15 A from its native structure) . Still, 
in agreement with hydrogen exchange studies (po|), some 
protection is expected across domains from our simula- 
tions and complete protection from exchange is expected 
only after the formation of the fully folded protein. A 
combination of fluorescence, CD mutagenic and new drug 
binding studies on DHFR indeed demonstrate that do- 
main 1 is largely folded with specific tertiary contacts 
formed and that this collection of intermediates is oblig- 
atory in the folding route (p^). 

B. Interleukin— 1/3 

Supported by some recent experiments, Heidary et al 
( p5| ) have proposed a kinetic mechanism for the folding 
of IL-1/3 that requires the presence of a well defined on- 
pathway intermediate species. The structural details of 
these species were determined from NMR and hydrogen 
exchange techniques (|2^j2^). We have compared these 
experimental data with our simulations for the IL-1/3 Go- 
like analogue ( Fig. |^ (b) and (d)). The folding picture 
emerging from these numerical studies differs substan- 
tially from that observed for DHFR (see panels (a) and 
(c) of Fig. |l|). An intermediate state is populated for Q 
around 0.55, followed by a rate limiting barrier (around 
Q = 0.7) after which the system proceeds to the well 
defined native state. 

Is the theoretical intermediate similar to the one ob- 
served experimentally? Using the same procedure em- 
ployed for the DHFR, a comparison between the average 



structure of the IL-1/3 intermediate ensemble and the one 
emerging from experimental studies is shown in Fig. IJ. 
These results indicate that the calculated intermediate 
has residues 40-105 (strands 4-8) folded into a native-like 
topology but with interactions between strands 5 and 8 
not fully completed. Experimental results confirm that 
strands 6-8 are well folded in the intermediate state and 
that strands 4-5 are partially formed. However results 
of experiments and theory differ in the region between 
residues 110-125 where hydrogen exchange shows early 
protection and theory predicts late contact formation. 
This region contains 4 aromatic groups PHE 112, PHE 
117, TYR 120, TRP 121 which may be sequestered from 
solvent due to clustering of these residues and from re- 
moval from unfavorable solvent interactions. This effect 
would not be fully accounted for our model, where all 
native interactions are considered as energetically equiv- 
alent and large stabilizing interactions are not differenti- 
ated. Thus, energetics may favor early formation of the 
structure corresponding to residues 105-125 while topol- 
ogy considerations favor the formation of strands 4-8. 

IV. CONCLUSIONS 

Theoretical and experimental studies of protein fold- 
ing at times appear to be at odds. Theoretical analysis of 
simple model systems oftentimes predict a large number 
of routes to the native protein whereas experimental work 
on larger systems indicates that folding proceeds through 
a limited number of intermediate species. Although in 
the eyes of some people these two descriptions are incon- 
sistent with each other, this is clearly not true. The large 
number of routes may or may not lead to the production 
of on-route kinetic intermediate ensembles depending on 
the result of the competition between configurational en- 
tropy and the effective folding energy. In this study, we 
show that productive intermediate species are produced 
by using simplified protein models, with funnel-like land- 
scapes, based on purely topological considerations and 
the results are in good agreement with the available ex- 
perimental data. The fact that these simplified minimally 
frustrated models for DHFR and IL-1/3 can predict the 
overall features of the folding intermediates and transi- 
tion states experimentally measured for these two pro- 
teins, with completely different folding mechanisms and 
functions, support our general picture that real proteins 
have a substantially reduced level of energetic frustra- 
tion and a large component of the observed heterogene- 
ity during the folding event is topologically determined. 
Such observations lead us to propose that the success 
in designing sequences that fold to a particular shape is 
constrained by topological effects. What is more chal- 
lenging are the consequences of this conclusion — are 
these topological constraints something that only have 
to be tolerated during the folding event or are they actu- 
ally used by biology towards helping function? Here we 
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speculate only in the context of these two examples, but 
this question really should be addressed more generally 
in the future. 
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FIGURE CAPTIONS 

Fig. 1 (a) RMS distances between the DHFR native 
structure and several computationally determined struc- 
tures at different values of the reaction coordinate Q for 
an unfolding simulation at a temperature slightly above 
the folding temperature (T = l.OlTy) and (c) free en- 
ergy F{Q) of the DHFR Go-like model as a function of 
Q around the folding temperature. The folding temper- 
ature Tf is estimated as the temperature where a sharp 
peak appears in the specific heat plotted as a function of 
the temperature (data not shown). Both temperatures 
and free energies are presented in units of Tf. Notice 
that the thermal fluctuations around the lowest energy 
state (i.e. Q = 1, by construction of the model) account 
for motions around the free energy minimum. There- 
fore, the folded state ensemble has a minimum close to 
Q — Q ^ 0.9, but not exactly at Q = 1. Indeed at 
Q — 1 the structure would be frozen in the native con- 
figuration. A similar remark applies for the IL-1/? free 
energy profile shown in panel (c). The energy of a con- 
figuration, as quantified by the color scale on the top of 
the figure, is here defined as the bare value of the effective 
potential function in that configuration (i.e. no configu- 
rational entropy is accounted in the energy) . Differences 
between energy and free energy (at finite temperature) 
are due to the configurational entropy contribution to the 
free energy. In panel (c) a main intermediate ensemble 
Ihf emerges in the folding process as a local minimum 
at Q around 0.4 after the overcoming of the first barrier. 
Indeed this local minimum corresponds to a populated 
region in panel (a) (after the scarcely populated barrier 
around Q = 0.3) with energy significantly lower than in 
the unfolded state. This main intermediate then evolves 
toward a set of structures close to the native states (lo- 
cated between Q = 0.7 and Q — 0.8) that eventually 
interconvert into the fully folded state. A transient set 
of structures, close to the native state, is also apparent 
from Fig. |^ (b) . The folding scheme resulting from these 
simulations is consistent with the sketch of Fig. ^ (d), 
proposed from the experimental data ( [l9| , |2^ ). 

(b) RMS distances between the native structure and 
several computationally determined structures at differ- 
ent values of the reaction coordinate Q for a folding sim- 
ulation of the Go-like model of IL-1/?. The simulation 
is performed at a temperature near to the folding tem- 
perature (T = O.gOT/). (d) Free energy F{Q) as a 
function of Q around to the folding temperature. The 
folding temperature is estimated from the sharp peak in 
the specific heat curve as a function of the temperature 
(data not shown). An intermediate ensemble is popu- 
lated during the folding event and it is identified by the 
broad local minimum in the free energy profile (around 
Q — 0.55), and the corresponding populated region in 
panel (b) (with energy significantly lower than in the 
unfolded state). These results are consistent with the 
kinetic mechanism for the folding of TL-1[3 proposed by 
Heidary et al. ([25|). A set of structures close to the na- 
tive conformation is transiently populated for Q between 



0.75 and 0.8 (see panel (b) and the corresponding "flat" 
region in the free energy panel (d) ). This fact could be 
interpreted as the presence of an additional intermedi- 
ate state close to the native state. Experimentally the 
possibility that another partially unfolded form could be 
populated during the folding process is currently under 
investigation. Several constant temperature simulations 
(both folding and unfolding simulations) of the two pro- 
tein models were made and combined to generate the free 
energy plots. 



Fig. 2 The probability Q,j{Q) of the native DHFR 
contacts to be formed, as resulting from the simulations 
at different stages of the folding process: (a) at an early 
stage {Q = 0.1 ± 0.05), (b) at the main intermediate - 
located in the interval Q = 0.4 ± 0.05 (see panels (a) 
and (c) of Fig. ^ and (c) at a late stage of the fold- 
ing process {Q = 0.7 ± 0.05). In an topologically and 
energetically perfectly smooth funnel-like energy land- 
scape, at any value Q during the folding, any contact 
(i, j) should have a probability Qij (Q) to be formed equal 
to Q (||). By computing Qij{Q) for each contact over 
different windows of the reaction coordinate Q, we can 
quantify the deviations from this smooth funnel behavior 
and locate the early and late contacts along the folding 
process. It is worth noticing that any deviation from the 
"perfectly smooth" behavior is mainly due to topological 
constraints, since energetic frustration has been mostly 
removed from the system. Different colors in the contact 
maps indicate different probability values from to 1, as 
quantified by the color scale at the top. The preference 
to form more local structure than non-local in the almost 
unfolded state (b) is due to the smaller conformational 
entropy loss by forming local contacts than by pinching 
off longer loops (^^ . The most interesting result is that 
domain 1, identified by the interactions among strands 
2-5 and helices 2-3, is substantially formed at the inter- 
mediate Ihf (probabilities for individual contacts grater 
than 0.7), while the formation of domain 2 (i.e. interac- 
tions among strand 1, strands 6-8, helix 1 and helix 4) 
is highly unfolded (contact probabilities between and 
0.4). Helix 1 and helix 4 are largely formed, but their 
interactions with the remainder of the proteins are loose 
(probabilities less than 0.4). Overall, this description of 
the structure of the main intermediate Ihf ~ domain 1 
almost formed and domain 2 largely unformed - is in 
agreement with the structure of Ihf experimentally ob- 
served. Moreover, the latest events in the folding process 
(panel (c)) appear to be the formation of interactions 
between strands 7-8 and the remainder of the protein. 
This again has been experimentally determined. Panel 
(d) illustrates the regions of the native structure that 
simulations and experiments agree to indicate as formed 
at the intermediate Ihf- 
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Fig. 3 The RMS distances between the regions of 
the DHFR structure identified as (a) domain 1 and (c) 
domain 2 and their corresponding native configurations 
are plotted versus the reaction coordinate Q. Domain 
1 collapses to a structure close to its native conforma- 
tion (RMS distance less than 5 A) in the early stages of 
folding, leading to the formation of the main intermediate 
Ihf (located at Q around 0.4) whereas domain ^remains 
largely unfolded (RMS distance larger than 15 A). In the 
interval of Q from 0.6 to 0.8 there are several possible 
structures. Consistently with the multi-channel folding 
model proposed from experimental evidences (19), from 
panels (a) and (c) one can propound several possible ways 
to proceed from Ihf to the folded state. In panel (b) 
the fraction of native contacts formed, Q, is plotted ver- 
sus the simulation time for a region of our simulations 
where the transition from folded to unfolded state is ob- 
served (at a simulation temperature slightly higher than 
the folding temperature, T — l.QlTj). A set of structures 
close to the native state {Q around 0.7) is transiently pop- 
ulated. Different colors represent different energies of a 
configuration (quantified by the top energy scale), as well 
as for panels (a) and (c). (d) The kinetics mechanism for 
the folding of DHFR, proposed on the basis of experimen- 
tal results (p^|-p^. Experimentally a first step of folding 
is detected as a very rapid collapse of the unfolded form 
to the burst-phase intermediate {I bp) which has a sig- 
nificant content of secondary structure. The folding state 
is reached through four different channel, involving the 
formations of the main intermediate Ihf- Ihf is repre- 
sented as a set of structures h — 14 structurally similar to 
each other but proceeding towards the native state with 
a different rate. These intermediate structures evolve to 
the native forms iVi — N4 via slow-folding reactions. 



(d) . The small difference between simulations and exper- 
iments (contact formation in the region between residues 
110-125) may be due to energetic considerations that are 
not taken into account in the model, as discussed in the 
text. In agreement with experimental results, the forma- 
tion of contacts between N and C terminus is not accom- 
plished until the late stage of folding - these contacts are 
still unformed for Q = (0.6 — 0.7) (data not shown). 



Fig. 4 Probability of the contact formation for the 
native contacts, as obtained during a typical folding sim- 
ulation (data shown are obtained at T = 0.99T/) of the 
IL-1/3 at different stages of the folding: (a) in a range 
of Q between 0.3 and 0.4 that corresponds to the early 
stage of folding leading to the formation of the inter- 
mediate ensemble; and (b) at the intermediate {Q be- 
tween 0.45 and 0.55). At the intermediate the interac- 
tions involving strands from 4 to 8 are almost completely 
formed; interactions among strands 1-3 are likely formed 
but interactions between them and the rest of the protein 
are loose. Contacts involving strands 9-12 appear weak- 
ened and the interactions between N (residues 1-40) and 
C terminus (residues 110-153) are completely unformed. 
Experimental results confirm that strands 6-8 are well 
folded in the intermediate state and that strands 4-5 are 
partially formed. Panels (c) and (d) show the regions 
of the IL-1/3 native structure formed at the intermedi- 
ate, as resulting from simulations (c) and experiments 
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FIG. 2. 




FIG. 3. 
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FIG. 4. 



11 



